import data from 86 files .nc4¶
- date: 2019-01-01
- variables:
- latitudes
- longitudes
- CO2_std
- CO2_sdev
- dimensions:
- time
- lat
- lon
- attributes:
- year
- month
- day
- time
- variables:
Normailize data ???¶
In [ ]:
# import libs
import os
import netCDF4 as nc4
import pandas as pd
import geopandas as gpd
from shapely.geometry import Point
import numpy as np
import itertools
In [ ]:
# Ruta al directorio que contiene los archivos netCDF
directorio_archivos = 'descargas'
# Lista de archivos netCDF en el directorio
archivos_netCDF = [os.path.join(directorio_archivos, nombre_archivo) for nombre_archivo in os.listdir(directorio_archivos) if nombre_archivo.endswith('.nc4')]
# create a dict de dimensiones [91,144]
row = np.arange(0,91)
col = np.arange(0,144)
dict_dim = dict()
def media_armonica(valores):
inversos = [1 / x for x in valores if x != 0] # Calcula los inversos de los valores no iguales a 0
if len(inversos) == 0:
return 0 # Si no hay valores válidos, devuelve 0
return len(inversos) / sum(inversos)
for i,j in itertools.product(row,col):
dict_dim[(i,j)] = {"lat": [], "lon": [], "co2_std": [], "co2_sdev": []}
# dict_dim[(i,j,"lon")] = []
# dict_dim[(i,j,"co2_std")] = []
# dict_dim[(i,j,"co2_sdev")] = []
for i in range(len(archivos_netCDF)):
# Abre el archivo netCDF
data = nc4.Dataset(archivos_netCDF[i],"r")
latitude = data.variables['Latitude'][:]
longitude = data.variables['Longitude'][:]
co2_std = data.variables['mole_fraction_of_carbon_dioxide_in_free_troposphere'][:]
co2_sdev = data.variables['mole_fraction_of_carbon_dioxide_in_free_troposphere_sdev'][:]
# print(latitude)
# print(longitude)
# print(co2_std)
# print(co2_sdev)
lon, lat = np.meshgrid(longitude, latitude)
# Aplana las mallas de coordenadas y la matriz de co2_free_troposphere a vectores unidimensionales
lon_flat = lon #.flatten()
lat_flat = lat #.flatten()
# create data.fream
df_lat = pd.DataFrame(lat_flat).fillna(0)
df_lon = pd.DataFrame(lon_flat).fillna(0)
df_co2_std = pd.DataFrame(co2_std).fillna(0)
df_co2_sdev = pd.DataFrame(co2_sdev).fillna(0)
for i,j in itertools.product(row,col):
aux_lat = df_lat.iloc[i,j]
dict_dim[(i,j)]["lat"].append(aux_lat)
dict_dim[(i,j)]["lon"].append(df_lon.iloc[i,j])
dict_dim[(i,j)]["co2_std"].append(df_co2_std.iloc[i,j])
dict_dim[(i,j)]["co2_sdev"].append(df_co2_sdev.iloc[i,j])
# Crear un data frame de los puntos unicos lat, lon para la dimension i,j
# print(dict_dim[(0,0)]["lat"])
columns = ["co2_std", "co2_sdev", "lat", "lon"]
df_dim = pd.DataFrame(columns=columns)
glat, glon = [], []
for i,j in itertools.product(row,col):
# print(dict_dim[(i,j)]["lat"])
# print(dict_dim[(i,j)]["lon"])
# print(dict_dim[(i,j)]["co2_std"])
# print(dict_dim[(i,j)]["co2_sdev"])
# print("==================================")
# print("==================================")
# print("==================================")
# print("=============================
x = media_armonica(dict_dim[(i,j)]["co2_sdev"])
y = media_armonica(dict_dim[(i,j)]["co2_std"])
#print(pd.Series(dict_dim[(i,j)]["lat"]).unique()[0], pd.Series(dict_dim[(i,j)]["lon"]).unique()[0])
glat.append(pd.Series(dict_dim[(i,j)]["lat"]).unique()[0]), glon.append(pd.Series(dict_dim[(i,j)]["lon"]).unique()[0])
min_lat_lon = min([len(pd.Series(dict_dim[(i,j)]["lat"]).unique()),len(pd.Series(dict_dim[(i,j)]["lon"]).unique())])
if min_lat_lon == 1:
df_dim.loc[len(glat)-1] = [x, y, pd.Series(dict_dim[(i,j)]["lat"]).unique()[0], pd.Series(dict_dim[(i,j)]["lon"]).unique()[0]]
# ({"co2_std": x, "co2_sdev": y, "lat": [pd.Series(dict_dim[(i,j)]["lat"]).unique()[0]], "lon": [pd.Series(dict_dim[(i,j)]["lon"]).unique()[0]]})
else:
df_dim.loc[len(glat)-1] = [x, y, pd.Series(dict_dim[(i,j)]["lat"]).unique()[0], pd.Series(dict_dim[(i,j)]["lon"]).unique()[0]]
print("no es unico".upper().ljust(30,"="))
# print()
# value.append(df_lat.iloc[i,j])
# value.append(df_lon.iloc[i,j])
In [ ]:
# normalize df_dim
In [ ]:
#create a geodataframe with geopandas and shapely: lat, lon -> glat, glon
#from matplotlib.patches import Polygon
#geometry = [gdp.Polygon(xy) for xy in zip(glon, glat)]
geometry = [Point(xy) for xy in zip(glon, glat)]
crs = 'EPSG:4326'
gdf = gpd.GeoDataFrame(df_dim, crs=crs, geometry=geometry)
In [ ]:
min(df_dim["co2_std"])
float(max(df_dim["co2_std"]))
df_dim.query("co2_std >= 0.000005")
Out[ ]:
| co2_std | co2_sdev | lat | lon | |
|---|---|---|---|---|
| 269 | 0.000005 | 0.000396 | 88.0 | 132.5 |
In [ ]:
#describe the data
gdf.describe()
Out[ ]:
| co2_std | co2_sdev | lat | lon | |
|---|---|---|---|---|
| count | 13104.000000 | 13104.000000 | 13104.000000 | 13104.000000 |
| mean | 0.000002 | 0.000330 | -0.005495 | -1.250000 |
| std | 0.000001 | 0.000148 | 52.528319 | 103.924508 |
| min | 0.000000 | 0.000000 | -90.000000 | -180.000000 |
| 25% | 0.000002 | 0.000394 | -46.000000 | -90.625000 |
| 50% | 0.000003 | 0.000395 | 0.000000 | -1.250000 |
| 75% | 0.000003 | 0.000396 | 46.000000 | 88.125000 |
| max | 0.000005 | 0.000419 | 89.500000 | 177.500000 |
In [ ]:
gdf.info()
<class 'geopandas.geodataframe.GeoDataFrame'> Index: 13104 entries, 0 to 13103 Data columns (total 5 columns): # Column Non-Null Count Dtype --- ------ -------------- ----- 0 co2_std 13104 non-null float64 1 co2_sdev 13104 non-null float64 2 lat 13104 non-null float64 3 lon 13104 non-null float64 4 geometry 13104 non-null geometry dtypes: float64(4), geometry(1) memory usage: 614.2 KB
In [ ]:
import plotly.express as px
fig = px.histogram(df_dim, x="co2_std", nbins=80, title="Histograma de co2_std")
fig.show()
In [ ]:
import plotly.express as px
fig = px.histogram(df_dim, x="co2_sdev", nbins=80, title="Histograma de co2_sdev")
fig.show()
Basado en el segso observado en los graficos, optamos por eliminar valores donde co2_std -> 0¶
In [ ]:
gdf2 = gdf.query("co2_std >= 0.000000001") #.plot(column='co2_std', cmap='OrRd', figsize=(15, 15), legend=True, scheme='quantiles')
In [ ]:
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib.colors as mcolors
import matplotlib.cm as cm
In [ ]:
gdf2["co2_std"]
Out[ ]:
8 0.000003
91 0.000003
119 0.000004
144 0.000002
145 0.000002
...
10939 0.000002
10940 0.000003
10941 0.000002
10942 0.000003
10943 0.000003
Name: co2_std, Length: 10803, dtype: float64
In [ ]:
# Crea una figura de Matplotlib
fig, ax = plt.subplots(figsize=(12, 8))
# Dibuja el mapa del mundo de GeoPandas como fondo
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.boundary.plot(ax=ax, linewidth=1, color='black')
# Utiliza Seaborn para dibujar los puntos y personalizar la paleta de colores
sns.scatterplot(ax=ax,x= 'lon', y= 'lat', data = gdf2, hue = 'co2_std', palette='coolwarm')
# Agrega opacidad a los puntos para diferenciar la concentración
alpha = 0.3
ax.collections[0].set_alpha(alpha)
# Personaliza el aspecto del mapa
ax.set_title('Distribución de CO2 en la Troposfera Libre con Opacidad')
plt.axis('off') # Desactiva los ejes
sm = cm.ScalarMappable(cmap='coolwarm')
sm.set_norm(mcolors.Normalize(vmin=min(gdf2["co2_std"]), vmax=max(gdf2["co2_std"])))
# Agrega una barra de colores (colorbar)
cbar = plt.colorbar(sm, ax=ax)
cbar.set_label('Concentración de CO2 (co2_std)')
# Muestra el mapa con los datos
plt.show()
/tmp/ipykernel_1443/3664182266.py:5: FutureWarning: The geopandas.dataset module is deprecated and will be removed in GeoPandas 1.0. You can get the original 'naturalearth_lowres' data from https://www.naturalearthdata.com/downloads/110m-cultural-vectors/.
In [ ]:
# Crea una figura de Matplotlib
fig, ax = plt.subplots(figsize=(12, 8))
# Dibuja el mapa del mundo de GeoPandas como fondo
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.boundary.plot(ax=ax, linewidth=1, color='black')
# Utiliza Seaborn para dibujar los puntos y personalizar la paleta de colores
sns.scatterplot(ax=ax,x= 'lon', y= 'lat', data = gdf2, hue = 'co2_sdev', palette='viridis')
# Agrega opacidad a los puntos para diferenciar la concentración
alpha = 0.3
ax.collections[0].set_alpha(alpha)
# Personaliza el aspecto del mapa
ax.set_title('Distribución de CO2_sdev en la Troposfera Libre con Opacidad')
plt.axis('off') # Desactiva los ejes
sm = cm.ScalarMappable(cmap='viridis')
sm.set_norm(mcolors.Normalize(vmin=min(gdf2["co2_sdev"]), vmax=max(gdf2["co2_sdev"])))
# Agrega una barra de colores (colorbar)
cbar = plt.colorbar(sm, ax=ax)
cbar.set_label('Concentración de CO2 (co2_sdev)')
# Muestra el mapa con los datos
plt.show()
/tmp/ipykernel_1443/1589977332.py:5: FutureWarning: The geopandas.dataset module is deprecated and will be removed in GeoPandas 1.0. You can get the original 'naturalearth_lowres' data from https://www.naturalearthdata.com/downloads/110m-cultural-vectors/.